EEG analysis

imports

library(tidyr)
library(pspline)
library(purrr)
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(reshape2)
library(TTR)
require(smooth)
require(greybox)
require(Mcomp)
library(RColorBrewer)
library(forecast)

add_loess <- function(df){
  loess_df <- data.frame(df$timestamps)
  for(i in names(df)){
    if(grepl('absolute', i) & !grepl('loess', i)){
      name <- paste("loess_", i, sep = "")
      df[, name] <- loess(df[,i] ~ df$timestamps, data = df, span=0.65)$fitted
    }
  }
  return(df)
}

add_average <- function(df, wave){
  name <- paste("avg_", wave, sep = "")
  temp_df <- df[,names(df) %in% colnames(df)[grepl(wave,colnames(df)) & !grepl("loess",colnames(df)) & !grepl("1",colnames(df)) & !grepl("4",colnames(df))]]
  temp_avg <- rowMeans(temp_df, dims = 1)
  df[, name] <- temp_avg
  return(df)
}

load data

inexp_meditation_files = sort(list.files("ffted/",pattern="^0_ffted_med"))
inexp_reference_files = sort(list.files("ffted/",pattern="^0_ffted_ref"))
exp_meditation_files = sort(list.files("ffted/",pattern="^1_ffted_med"))
exp_reference_files = sort(list.files("ffted/",pattern="^1_ffted_ref"))
exp_meditation = list()
for(i in 1:length(exp_meditation_files)) {
  file = exp_meditation_files[i]
  exp_meditation[[i]] <- read.csv(paste("ffted/", file, sep=""))
  exp_meditation[[i]] <- exp_meditation[[i]][rowSums(exp_meditation[[i]] == "-Inf") == 0, , drop = FALSE]
  exp_meditation[[i]] <- add_average(exp_meditation[[i]], "alpha")
  exp_meditation[[i]] <- add_average(exp_meditation[[i]], "beta")
  exp_meditation[[i]] <- add_average(exp_meditation[[i]], "gamma")
  exp_meditation[[i]] <- add_average(exp_meditation[[i]], "delta")
  exp_meditation[[i]] <- add_average(exp_meditation[[i]], "theta")
  exp_meditation[[i]] <- add_loess(exp_meditation[[i]])
}
inexp_meditation = list()
for(i in 1:length(inexp_meditation_files)) {
  file = inexp_meditation_files[i]
  inexp_meditation[[i]] <-  read.csv(paste("ffted/", file, sep=""))
  inexp_meditation[[i]] <- inexp_meditation[[i]][rowSums(inexp_meditation[[i]] == "-Inf") == 0, , drop = FALSE]
  inexp_meditation[[i]] <- add_average(inexp_meditation[[i]], "alpha")
  inexp_meditation[[i]] <- add_average(inexp_meditation[[i]], "beta")
  inexp_meditation[[i]] <- add_average(inexp_meditation[[i]], "gamma")
  inexp_meditation[[i]] <- add_average(inexp_meditation[[i]], "delta")
  inexp_meditation[[i]] <- add_average(inexp_meditation[[i]], "theta")
  inexp_meditation[[i]] <- add_loess(inexp_meditation[[i]])
}
exp_reference = list()
for(i in 1:length(exp_reference_files)) {
  file = exp_reference_files[i]
  exp_reference[[i]] <-  read.csv(paste("ffted/", file, sep=""))
  exp_reference[[i]] <- exp_reference[[i]][rowSums(exp_reference[[i]] == "-Inf") == 0, , drop = FALSE]
  exp_reference[[i]] <- add_average(exp_reference[[i]], "alpha")
  exp_reference[[i]] <- add_average(exp_reference[[i]], "beta")
  exp_reference[[i]] <- add_average(exp_reference[[i]], "gamma")
  exp_reference[[i]] <- add_average(exp_reference[[i]], "delta")
  exp_reference[[i]] <- add_average(exp_reference[[i]], "theta")
  exp_reference[[i]] <- add_loess(exp_reference[[i]])
}
inexp_reference = list()
for(i in 1:length(inexp_reference_files)) {
  file = inexp_reference_files[i]
  inexp_reference[[i]] <-  read.csv(paste("ffted/", file, sep=""))
  inexp_reference[[i]] <- inexp_reference[[i]][rowSums(inexp_reference[[i]] == "-Inf") == 0, , drop = FALSE]
  inexp_reference[[i]] <- add_average(inexp_reference[[i]], "alpha")
  inexp_reference[[i]] <- add_average(inexp_reference[[i]], "beta")
  inexp_reference[[i]] <- add_average(inexp_reference[[i]], "gamma")
  inexp_reference[[i]] <- add_average(inexp_reference[[i]], "delta")
  inexp_reference[[i]] <- add_average(inexp_reference[[i]], "theta")
  inexp_reference[[i]] <- add_loess(inexp_reference[[i]])
}

convert time from absolute to relative

convert_time <- function(timestamps){
  time <- timestamps - min(timestamps)
  return(time)
}

subsetting data

subset_and_melt <- function(med, exp, wave, electrodes, melt, melt_all){
  if(med){
    if(exp){
      df <- exp_meditation
    }
    else{
      df <- inexp_meditation
    }
  }
  else{
    if(exp){
      df <- exp_reference
    }
    else{
      df <- inexp_reference
    }
  }
  data = list()
  for(i in 1:length(df)) {
  data[[i]] <- df[[i]][, grepl(paste('(', wave, ')|(', electrodes, ')|(timestamp)', sep = ""), names(df[[i]]))]
  data[[i]]$timestamps <- convert_time(data[[i]]$timestamps)
  if(melt){
    data[[i]] <- melt(data[[i]], id.vars = 'timestamps', variable.name = 'waves')
    data[[i]]$waves <- as.factor(data[[i]]$waves)
  }
  }
  if(melt_all){
    data <- melt(data, id.vars = c('timestamps', 'waves', 'value'))
    data$L1 <- as.factor(data$L1)
  }
  
  return(data)
}

testing

temp_2 <- subset_and_melt(FALSE, TRUE, "beta", "None", TRUE, TRUE)
#ggplot(temp_2[[1]], aes(timestamps,value)) + geom_line(aes(colour = waves))

things to think about: * extreme points, spikes, etc * different smoothing rate * difference between meditative and regular state * difference between experienced meditators and newbies * how to determine the levels and scale it for others * what happens when meditation is stable * explore the coherence between alpha and theta waves

explore experienced/inexperienced distributions of alpha wave

alpha_exp_med = subset_and_melt(TRUE, TRUE, 'alpha', 'ALL', TRUE, TRUE)
alpha_inexp_med = subset_and_melt(TRUE, FALSE, 'alpha', 'ALL', TRUE, TRUE)
alpha_exp_med$exp <- "experienced"
alpha_exp_med$exp <- as.factor(alpha_exp_med$exp)
alpha_inexp_med$exp <- "inexperienced"
alpha_inexp_med$exp <- as.factor(alpha_inexp_med$exp)
temp <- rbind(alpha_exp_med, alpha_inexp_med)
ggplot(temp, aes(value, fill = exp)) + 
  geom_histogram(alpha = 0.5, bins = 100, position = 'identity')

explore experienced/inexperienced distributions of theta wave

theta_exp_med = subset_and_melt(TRUE, TRUE, 'theta', 'ALL', TRUE, TRUE)
theta_inexp_med = subset_and_melt(TRUE, FALSE, 'theta', 'ALL', TRUE, TRUE)
theta_exp_med$exp <- "experienced"
theta_exp_med$exp <- as.factor(theta_exp_med$exp)
theta_inexp_med$exp <- "inexperienced"
theta_inexp_med$exp <- as.factor(theta_inexp_med$exp)
temp <- rbind(theta_exp_med, theta_inexp_med)
ggplot(temp, aes(value, fill = exp)) + 
  geom_histogram(alpha = 0.5, bins = 100, position = 'identity')

explore experienced/inexperienced distributions of beta wave

beta_exp_med = subset_and_melt(TRUE, TRUE, 'beta', 'ALL', TRUE, TRUE)
beta_inexp_med = subset_and_melt(TRUE, FALSE, 'beta', 'ALL', TRUE, TRUE)
beta_exp_med$exp <- "experienced"
beta_exp_med$exp <- as.factor(beta_exp_med$exp)
beta_inexp_med$exp <- "inexperienced"
beta_inexp_med$exp <- as.factor(beta_inexp_med$exp)
temp <- rbind(beta_exp_med, beta_inexp_med)
ggplot(temp, aes(value, fill = exp)) + 
  geom_histogram(alpha = 0.5, bins = 100, position = 'identity')

explore experienced/inexperienced distributions of gamma wave

gamma_exp_med = subset_and_melt(TRUE, TRUE, 'gamma', 'ALL', TRUE, TRUE)
gamma_inexp_med = subset_and_melt(TRUE, FALSE, 'gamma', 'ALL', TRUE, TRUE)
gamma_exp_med$exp <- "experienced"
gamma_exp_med$exp <- as.factor(gamma_exp_med$exp)
gamma_inexp_med$exp <- "inexperienced"
gamma_inexp_med$exp <- as.factor(gamma_inexp_med$exp)
temp <- rbind(gamma_exp_med, gamma_inexp_med)
ggplot(temp, aes(value, fill = exp)) + 
  geom_histogram(alpha = 0.5, bins = 100, position = 'identity')

explore experienced/inexperienced distributions of delta wave

t.test(na.omit(alpha_exp_med$value), na.omit(alpha_inexp_med$value), var.equal = TRUE)

    Two Sample t-test

data:  na.omit(alpha_exp_med$value) and na.omit(alpha_inexp_med$value)
t = 101.66, df = 685930, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2013323 0.2092484
sample estimates:
mean of x mean of y 
0.9737287 0.7684384 

t.test(na.omit(alpha_exp_med$value), na.omit(alpha_inexp_med$value), var.equal = TRUE)

    Two Sample t-test

data:  na.omit(alpha_exp_med$value) and na.omit(alpha_inexp_med$value)
t = 161.49, df = 1543400, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2040037 0.2090164
sample estimates:
mean of x mean of y 
0.9426494 0.7361393 

t_test <- function(wave, med){
  wave_exp <- subset_and_melt(med, TRUE, wave, "All", TRUE, FALSE)
  medians_exp <- list()
  for(i in 1:length(wave_exp)){
    medians_exp[[i]] <- median(wave_exp[[i]]$value)
  }
  wave_inexp <- subset_and_melt(med, FALSE, wave, "All", TRUE, FALSE)
  medians_inexp <- list()
  for(i in 1:length(wave_inexp)){
    medians_inexp[[i]] <- median(wave_inexp[[i]]$value)
  }
  
  medians_exp <- unlist(medians_exp, use.names=FALSE)
  medians_inexp <- unlist(medians_inexp, use.names=FALSE)
  return(print(t.test(medians_exp, medians_inexp)))
}

a <- t_test("beta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -0.18003, df = 7.7255, p-value = 0.8618
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.4753273  0.4068778
sample estimates:
mean of x mean of y 
0.8503494 0.8845742 

a <- t_test("gamma", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -0.43366, df = 7.0362, p-value = 0.6775
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.5995961  0.4135896
sample estimates:
mean of x mean of y 
0.4898110 0.5828143 

a <- t_test("delta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.76505, df = 9.4431, p-value = 0.4629
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2353533  0.4785046
sample estimates:
mean of x mean of y 
0.3027340 0.1811583 

a <- t_test("theta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.57218, df = 12.632, p-value = 0.5772
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2152788  0.3697712
sample estimates:
mean of x mean of y 
0.7990211 0.7217749 

a <- t_test("theta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.51721, df = 12.545, p-value = 0.614
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2674452  0.4349945
sample estimates:
mean of x mean of y 
0.7675215 0.6837468 

There is no significant statistical difference in medians and means between experienced and inexperienced people

t_test_min <- function(wave, med){
  wave_exp <- subset_and_melt(med, TRUE, wave, "All", TRUE, FALSE)
  medians_exp <- list()
  for(i in 1:length(wave_exp)){
    medians_exp[[i]] <- min(wave_exp[[i]]$value)
  }
  wave_inexp <- subset_and_melt(med, FALSE, wave, "All", TRUE, FALSE)
  medians_inexp <- list()
  for(i in 1:length(wave_inexp)){
    medians_inexp[[i]] <- min(wave_inexp[[i]]$value)
  }
  
  medians_exp <- unlist(medians_exp, use.names=FALSE)
  medians_inexp <- unlist(medians_inexp, use.names=FALSE)
  return(print(t.test(medians_exp, medians_inexp)))
}

a <- t_test_min("alpha", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.40198, df = 9.0746, p-value = 0.697
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.8983212  1.2871653
sample estimates:
mean of x mean of y 
-1.197179 -1.391601 

a <- t_test_min("theta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.47281, df = 9.3345, p-value = 0.6472
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.8883717  1.3610984
sample estimates:
mean of x mean of y 
-1.227769 -1.464133 

a <- t_test_min("gamma", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.1358, df = 7.1696, p-value = 0.8957
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.9793947  1.0993536
sample estimates:
 mean of x  mean of y 
-0.8495856 -0.9095651 

a <- t_test_min("beta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = 0.34185, df = 8.3868, p-value = 0.7409
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.007439  1.361425
sample estimates:
 mean of x  mean of y 
-0.4808075 -0.6578005 

a <- t_test_min("delta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -0.17516, df = 12.919, p-value = 0.8637
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.486542  1.263703
sample estimates:
mean of x mean of y 
-4.062885 -3.951465 

t_test_max <- function(wave, med){
  wave_exp <- subset_and_melt(med, TRUE, wave, "All", TRUE, FALSE)
  medians_exp <- list()
  for(i in 1:length(wave_exp)){
    medians_exp[[i]] <- max(wave_exp[[i]]$value)
  }
  wave_inexp <- subset_and_melt(med, FALSE, wave, "All", TRUE, FALSE)
  medians_inexp <- list()
  for(i in 1:length(wave_inexp)){
    medians_inexp[[i]] <- max(wave_inexp[[i]]$value)
  }
  
  medians_exp <- unlist(medians_exp, use.names=FALSE)
  medians_inexp <- unlist(medians_inexp, use.names=FALSE)
  return(print(t.test(medians_exp, medians_inexp)))
}

a <- t_test_max("alpha", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -1.0657, df = 11.867, p-value = 0.3078
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.2699213  0.4363986
sample estimates:
mean of x mean of y 
 3.694748  4.111509 

a <- t_test_max("theta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -0.56219, df = 10.179, p-value = 0.5862
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.9047506  0.5394850
sample estimates:
mean of x mean of y 
 4.309198  4.491831 

a <- t_test_max("gamma", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -1.7215, df = 9.7672, p-value = 0.1166
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.6313428  0.2118538
sample estimates:
mean of x mean of y 
 3.275454  3.985199 

a <- t_test_max("beta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -1.5388, df = 8.1208, p-value = 0.1619
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.4203471  0.2816698
sample estimates:
mean of x mean of y 
 3.641218  4.210557 

a <- t_test_max("delta", TRUE)

    Welch Two Sample t-test

data:  medians_exp and medians_inexp
t = -0.13081, df = 12.048, p-value = 0.8981
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.9272460  0.8221727
sample estimates:
mean of x mean of y 
 4.311494  4.364031 

#ggplot(temp, aes(value, fill = exp)) + 
#  geom_histogram(alpha = 0.5, bins = 100, position = 'identity')
data_temp <- subset_and_melt(TRUE, TRUE, 'alpha', "ALL", TRUE, TRUE)
#geom_line(aes(y=theta_absolute_2, x = timestamps), 
#       data = without_na, color=brewer.pal(4, "Blues")[3]) + #geom_smooth(aes(y=theta_absolute_2, x = timestamps), 
#       data = without_na, color=brewer.pal(4, "Blues")[3], span = 0.01) +

#temps <- data_temp[data_temp$L1 == 1 & data_temp$waves == "alpha_absolute_1", ]
# data_temp$L1 == 1,
#data_temp$waves == "alpha_absolute_1", 
#data_temp <- data_temp[data_temp$L1 == 1,]
to_draw <- data_temp[data_temp$L1 == 1 & data_temp$waves == "alpha_absolute_1", ]
#ggplot(aes(y = value, x = timestamps), data = data_temp, data_temp$waves == "alpha_absolute_1") +
#ggplot() + 
#  geom_line(aes(SMA(to_draw$value, n=15), x = to_draw$timestamps))
#ggplot(aes(y = value, x = timestamps), data = data_temp, data_temp$waves == "alpha_absolute_1") +
  #geom_line() +
#  geom_smooth(span = 0.001, level = 0.95, method = 'loess') 

temp <- es(to_draw$value, h=18, holdout=TRUE, silent=FALSE)

The provided data is not ts object. Only non-seasonal models are available.
Only additive models are allowed with non-positive data.

Forming the pool of models based on... ANN, AAN, Estimation progress:    100%... Done! 

#to_draw_part$sma <- sma(to_draw_part$value, n = 5, v = 0.9)$fitted
to_draw_part$loess <- loess(value ~ timestamps, data=to_draw_part, span=0.65)$fitted

Error in is.data.frame(data) : object 'to_draw_part' not found

ggplot() +
  geom_line(aes(y = value, x = timestamps), data = to_draw_part) +
  geom_smooth(aes(y = value, x = timestamps), data = to_draw_part, span = 1, n = 15, color = "blue") +
  geom_line(aes(y = loess, x = timestamps), data = to_draw_part, color = 'red')

ggplot() +
  geom_smooth(aes(y = value, x = timestamps), data = theta_exp_med[theta_exp_med$waves == 'loess_theta_absolute_2',], span = 1, n = 15, color = "blue") +
  geom_line(aes(y = value, x = timestamps), data = theta_exp_med[theta_exp_med$waves == 'loess_theta_absolute_2',], color = 'red', alpha = 0.2) +
  geom_line(aes(y = value, x = timestamps), data = theta_exp_med[theta_exp_med$waves == 'theta_absolute_2',], color = 'black', alpha = 0.2)

curr_data <- theta_exp_med[theta_exp_med$L1 == '1',]
#[theta_exp_med$waves == 'theta_absolute_2' & 
    
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'red', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1)

temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
loess_ordered <- temp$value
loess_ordered <- sort(loess_ordered)
min(loess_ordered)

[1] 0.3447938

barplot(loess_ordered)

temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

NA

theta_2, second person

curr_data <- theta_exp_med[theta_exp_med$L1 == '2',]
temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

curr_data <- theta_exp_med[theta_exp_med$L1 == '3',]
temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

curr_data <- theta_exp_med[theta_exp_med$L1 == '4',]
temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

curr_data <- theta_exp_med[theta_exp_med$L1 == '6',]
temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

curr_data <- theta_exp_med[theta_exp_med$L1 == '7',]
temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

curr_data <- theta_exp_med[theta_exp_med$L1 == '8',]
temp <- curr_data[curr_data$waves == 'loess_theta_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'theta_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_theta_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

curr_data <- alpha_exp_med[theta_exp_med$L1 == '1',]
temp <- curr_data[curr_data$waves == 'loess_alpha_absolute_2', ]
temp$state <- cut(temp$value, quantile(temp$value,(0:5)/5))
ggplot(aes(y = value, x = timestamps), data = curr_data) +
  geom_line(data = curr_data[curr_data$waves == 'alpha_absolute_2', ], color = 'grey', alpha = 0.7) +
  geom_line(data = curr_data[curr_data$waves == 'loess_alpha_absolute_2',], color = 'black', alpha = 1) +
  geom_line(aes(y = value, x = timestamps, color = state, size = 2),data = temp, alpha = 1) 

alpha_2_all_exp <- subset_and_melt(TRUE, TRUE, 'alpha', 'All', TRUE, TRUE)
ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_2',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 2, exp")

alpha_2_all_exp <- subset_and_melt(TRUE, TRUE, 'alpha', 'All', TRUE, TRUE)
ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_1',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 1, exp")

alpha_2_all_exp <- subset_and_melt(TRUE, TRUE, 'alpha', 'All', TRUE, TRUE)
ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_3',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 3, exp")

alpha_2_all_exp <- subset_and_melt(TRUE, TRUE, 'alpha', 'All', TRUE, TRUE)
ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_4',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 4, exp")

alpha_2_all_exp <- subset_and_melt(TRUE, FALSE, 'alpha', 'All', TRUE, TRUE)
ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_2',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 2, inexp")

ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_1',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 1, inexp")

ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_3',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 3, inexp")

ggplot(aes(x = timestamps, y = value), data = alpha_2_all_exp[alpha_2_all_exp$waves == 'loess_alpha_absolute_4',]) +
  geom_line(aes(color = L1)) +
  ggtitle("alpha 4, inexp")

temp <- subset_and_melt(TRUE, TRUE, 'theta', 'All', TRUE, TRUE)
ggplot(aes(x = timestamps, y = value), data = temp[temp$waves == 'loess_theta_absolute_1',]) +
  geom_line(aes(color = L1)) +
  ggtitle("theta 1, exp")

ggplot(aes(x = timestamps, y = value), data = temp[temp$waves == 'loess_theta_absolute_2',]) +
  geom_line(aes(color = L1)) +
  ggtitle("theta 2, exp")

ggplot(aes(x = timestamps, y = value), data = temp[temp$waves == 'loess_theta_absolute_3',]) +
  geom_line(aes(color = L1)) +
  ggtitle("theta 3, exp")

ggplot(aes(x = timestamps, y = value), data = temp[temp$waves == 'loess_theta_absolute_4',]) +
  geom_line(aes(color = L1)) +
  ggtitle("theta 4, exp")

temp <- subset_and_melt(TRUE, TRUE, 'alpha', 'theta', TRUE, TRUE)
ggplot() +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_theta_absolute_1',],  size = 1) +
  ggtitle("theta 1 and alpha 1, exp") +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_alpha_absolute_1',], linetype = 'dotted', size = 1) +
  scale_color_brewer(palette = "Spectral")

ggplot() +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_theta_absolute_2',],  size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_alpha_absolute_2',], linetype = 'dotted', size = 1) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("theta 2 and alpha 2, exp")

ggplot() +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_theta_absolute_3',],  size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_alpha_absolute_3',], linetype = 'dotted', size = 1) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("theta 3 and alpha 3, exp")

ggplot() +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_theta_absolute_4',],  size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = temp[temp$waves == 'loess_alpha_absolute_4',], linetype = 'dotted', size = 1) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("theta 4 and alpha 4, exp")

alpha_with_averages <- subset_and_melt(TRUE, TRUE, 'alpha', 'all', TRUE, TRUE)
ggplot() +
  geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'alpha_absolute_2',],  size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'loess_alpha_absolute_2',], linetype = 'dotted', size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], linetype = 'dotted', size = 1) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("alpha 2 all together, exp")

ggplot() +
  #geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == #'loess_alpha_absolute_2',], linetype = 'dotted', size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], size = 1, alpha = 0.5) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("alpha averages")

ggplot() +
  #geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == #'loess_alpha_absolute_2',], linetype = 'dotted', size = 1) +
  geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], size = 1, alpha = 0.2) +
  scale_color_brewer(palette = "Spectral") +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], size = 1) +
  ggtitle("alpha averages with smoothing")

ggplot() +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], size = 0.5) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("alpha averages with smoothing")

ggplot() +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], size = 0.5) +
    geom_line(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'loess_alpha_absolute_2',], size = 1, alpha = 1) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("alpha averages with smoothing VS alpha 2 (bold)")

theta_with_averages <- subset_and_melt(TRUE, TRUE, 'theta', 'all', TRUE, TRUE)

ggplot() +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = theta_with_averages[theta_with_averages$waves == 'avg_theta',], size = 0.5) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("theta averages with smoothing")

ggplot() +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = theta_with_averages[theta_with_averages$waves == 'avg_theta',], size = 0.5) +
    geom_line(aes(x = timestamps, y = value, color = L1), data = theta_with_averages[theta_with_averages$waves == 'loess_theta_absolute_2',], size = 1, alpha = 1) +
  scale_color_brewer(palette = "Spectral") +
  ggtitle("theta averages with smoothing VS theta 2 (bold)")

ggplot() +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = theta_with_averages[theta_with_averages$waves == 'avg_theta',], size = 0.5) +
  geom_smooth(aes(x = timestamps, y = value, color = L1), data = alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',], size = 1) +
  
  scale_color_brewer(palette = "Spectral") +
  ggtitle("theta averages VS alpha averages (bold)")

split_stages <- function(df, time_interval, representative_period, threshold){
  current_time <- 0
  current_stage <- 1
  df$stages <- rep(0,nrow(df))
  df$stages[df$timestamps < current_time + time_interval] <- current_stage #first 3 minutes
  while(TRUE){
    if(current_time + time_interval + 1 > max(df$timestamps)){
      break
    }
    repr_df <- df[(df$timestamps < (current_time + time_interval) & df$timestamps > (current_time + time_interval - representative_period)),]
    new_stage_sign <- next_stage(repr_df, threshold) #defining a new stage level
    if(new_stage_sign > 0){
      current_stage <- min(current_stage + 1, 5)
    }
    if(new_stage_sign < 0){
      current_stage <- max(current_stage - 1, 1)
    }
    #set a new stage in the df
    df[(df$timestamps < (current_time + 2 * time_interval) & df$timestamps > (current_time + time_interval)),]$stages <- current_stage
    current_time <- current_time + time_interval
  }
  df$stages <- as.factor(df$stages)
  return(df$stages)
}

next_stage <- function(df, threshold){
  fit <- auto.arima(df$value)
  forecast <- predict(fit, length(df$value))
  initial_data <- df[(df$value < quantile(df$value, 0.95) & df$value > quantile(df$value, 0.05)), ]$value
  forecast <- forecast$pred
  forecast <- forecast[forecast < quantile(forecast, 0.95) & forecast > quantile(forecast, 0.05)]
  #fit <- which(fit < quantile(fit, 0.95) & fit > quantile(fit, 0.05))
  diff <- sum(forecast) - sum(initial_data) 
  print(sum(forecast))
  print(sum(initial_data))
  print("-------")
  if(abs(diff) > threshold){
    if(diff > 0){
      return(1)
    }
    else{
      return(-1)
    }
  }
  return(0)
}

temp <- alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',]
temp <- temp[temp$L1 == '8',]
res <- split_stages(temp, 3 * 60, 30, 0.01)
print(unique(res))

[1] 1 2 3 4
Levels: 1 2 3 4

temp <- alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',]
temp <- temp[temp$L1 == '1',]
res <- split_stages(temp, 3 * 60, 30, 0.01)

add_stages <- function(df, time_interval, representative_period, threshold){
  df$stages <- rep(0,nrow(df))
  for(factor in levels(df$L1)){
    df[df$L1 == factor,]$stages <- split_stages( df[df$L1 == factor,], time_interval, representative_period, threshold)
  }
  return(df)
}

temp <- alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',]
temp <- add_stages(temp, 3 * 60, 30, 0.001)

alpha_with_averages <- subset_and_melt(TRUE, TRUE, 'alpha', 'all', TRUE, TRUE)
temp <- alpha_with_averages[alpha_with_averages$waves == 'avg_alpha',]
alpha_stages <- add_stages(temp, 3 * 60, 30, 5)

[1] -35.20688
[1] -34.29825
[1] "-------"
[1] -14.63786
[1] -15.28299
[1] "-------"
[1] 15.32103
[1] 10.06241
[1] "-------"
[1] 6.146795
[1] 82.95732
[1] "-------"
[1] 35.15567
[1] 36.35809
[1] "-------"
[1] 391.2204
[1] 390.8586
[1] "-------"
[1] 571.8857
[1] 536.5056
[1] "-------"
[1] 420.7765
[1] 422.2623
[1] "-------"
[1] 517.0316
[1] 498.9921
[1] "-------"
[1] 366.9362
[1] 416.9336
[1] "-------"
[1] 451.5398
[1] 440.67
[1] "-------"
[1] 468.0502
[1] 484.8621
[1] "-------"
[1] 510.2569
[1] 460.9515
[1] "-------"
[1] 48.10229
[1] 133.7512
[1] "-------"
[1] 98.76786
[1] 97.1298
[1] "-------"
[1] 84.62468
[1] 78.74971
[1] "-------"
[1] 86.0356
[1] 92.70745
[1] "-------"
[1] 5.766955
[1] 58.91794
[1] "-------"
[1] 71.68597
[1] 56.37134
[1] "-------"
[1] 436.6874
[1] 118.4694
[1] "-------"
[1] 81.68036
[1] 81.26119
[1] "-------"
[1] 64.96005
[1] 97.23356
[1] "-------"
[1] 81.21402
[1] 51.12887
[1] "-------"
[1] -2.428206
[1] 71.86342
[1] "-------"
[1] 49.71121
[1] 73.13329
[1] "-------"
[1] 84.41539
[1] 82.1399
[1] "-------"
[1] 67.44729
[1] 91.86051
[1] "-------"
[1] 65.09052
[1] 79.8176
[1] "-------"
[1] 122.6536
[1] 54.5804
[1] "-------"
[1] 57.36476
[1] 71.7172
[1] "-------"
[1] 0.01791189
[1] -1.269861
[1] "-------"
[1] -13.79072
[1] -0.02891286
[1] "-------"
[1] 9.225325
[1] 5.07441
[1] "-------"
[1] 5.871108
[1] -14.28683
[1] "-------"
[1] -0.01042182
[1] 4.490672
[1] "-------"
[1] -19.42542
[1] 35.55241
[1] "-------"
[1] -6.670456
[1] -0.07772321
[1] "-------"
[1] -7.039187e-05
[1] 4.346686
[1] "-------"
[1] 64.52856
[1] -4.243069
[1] "-------"
[1] -0.0002926155
[1] 8.918285
[1] "-------"
[1] 534.1726
[1] 532.2413
[1] "-------"
[1] 496.291
[1] 497.9683
[1] "-------"
[1] 473.2478
[1] 473.1004
[1] "-------"
[1] 194.3603
[1] 463.7927
[1] "-------"
[1] 475.0435
[1] 475.4863
[1] "-------"
[1] 459.2095
[1] 456.1062
[1] "-------"
[1] 465.478
[1] 451.3855
[1] "-------"
[1] 419.1453
[1] 455.8208
[1] "-------"
[1] 451.4211
[1] 450.9958
[1] "-------"
[1] 418.1923
[1] 444.9182
[1] "-------"

the_person <- alpha_stages[alpha_stages$L1 == '1',]
print(unique(the_person$stages))

[1] 1 2

ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '2',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '3',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '4',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '5',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '6',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '7',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")

the_person <- alpha_stages[alpha_stages$L1 == '8',]
ggplot() +
  geom_line(aes(y = value, x = timestamps, color = as.factor(stages)), data = the_person, alpha = 1) +
  geom_smooth(aes(y = value, x = timestamps), data = the_person, color = "black") +
  scale_x_continuous(breaks = seq(0, 1000, by = 180)) +
  scale_color_brewer(palette = "Blues")